!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine integrate (current)
!==================================================================================================================================|

   use parameters
   use block_define

   implicit none

   type(block), pointer, intent(inout) :: current

   integer(4) :: parent_id_p1, loop, i_loop, i, j
   real(8), dimension(nx_block, ny_block) :: ro_old, mx_old, my_old, mz_old, bx_old, by_old, bz_old, en_old, az_old, po_old
   real(8) :: cd_l
   type(block), pointer :: parent

!----------------------------------------------------------------------------------------------------------------------------------|
   if (current%nc .eq. 0) then

      cd_l = cd / current%lv

      loop = 1
      if (ns .eq. 1 .and. gravity_reset_flag .eq. 1) then
         loop = 3
         do j = 1, ny_block
         do i = 1, nx_block
            ro_old(i, j) = current%ro(i, j)
            mx_old(i, j) = current%mx(i, j)
            my_old(i, j) = current%my(i, j)
            mz_old(i, j) = current%mz(i, j)
            bx_old(i, j) = current%bx(i, j)
            by_old(i, j) = current%by(i, j)
            bz_old(i, j) = current%bz(i, j)
            en_old(i, j) = current%en(i, j)
            az_old(i, j) = current%az(i, j)
            po_old(i, j) = current%po(i, j)
         enddo
         enddo
      endif

      do i_loop = 1, loop
         select case(scheme)
         case (1)
            call tvd_mmc (current%ro, current%mx, current%my, current%mz,                                                          &
               current%bx, current%by, current%bz, current%en, current%az, current%po, current%gx, current%gy,                     &
               current%x, current%y,                                                                                               &
               nx_block, ny_block, ng, n_neighs, current%dx, current%dy, min_dt, max_et, gm, cd, cfl,                              &
               div_clean_flag, resistivity_flag,                                                                                   &
               current%fxro_bnd, current%fxmx_bnd, current%fxmy_bnd, current%fxmz_bnd,                                             &
               current%fxbx_bnd, current%fxby_bnd, current%fxbz_bnd, current%fxen_bnd, current%fxpo_bnd,                           &
               current%fyro_bnd, current%fymx_bnd, current%fymy_bnd, current%fymz_bnd,                                             &
               current%fybx_bnd, current%fyby_bnd, current%fybz_bnd, current%fyen_bnd, current%fypo_bnd)

         case (2)
            call tvd_lf (current%ro, current%mx, current%my, current%mz,                                                           &
               current%bx, current%by, current%bz, current%en, current%az, current%po, current%gx, current%gy,                     &
               current%x, current%y,                                                                                               &
               nx_block, ny_block, ng, n_neighs, current%dx, current%dy, min_dt, max_et, gm, cd, cfl,                              &
               div_clean_flag, resistivity_flag, riemann_solver_flag,                                                              &
               current%fxro_bnd, current%fxmx_bnd, current%fxmy_bnd, current%fxmz_bnd,                                             &
               current%fxbx_bnd, current%fxby_bnd, current%fxbz_bnd, current%fxen_bnd, current%fxpo_bnd,                           &
               current%fyro_bnd, current%fymx_bnd, current%fymy_bnd, current%fymz_bnd,                                             &
               current%fybx_bnd, current%fyby_bnd, current%fybz_bnd, current%fyen_bnd, current%fypo_bnd)

         case (3)
            call weno (current%ro, current%mx, current%my, current%mz,                                                             &
               current%bx, current%by, current%bz, current%en, current%az, current%po, current%gx, current%gy,                     &
               current%x, current%y,                                                                                               &
               nx_block, ny_block, ng, n_neighs, current%dx, current%dy, min_dt, max_et, gm, cd, cfl,                              &
               div_clean_flag, resistivity_flag, riemann_solver_flag,                                                              &
               current%fxro_bnd, current%fxmx_bnd, current%fxmy_bnd, current%fxmz_bnd,                                             &
               current%fxbx_bnd, current%fxby_bnd, current%fxbz_bnd, current%fxen_bnd, current%fxpo_bnd,                           &
               current%fyro_bnd, current%fymx_bnd, current%fymy_bnd, current%fymz_bnd,                                             &
               current%fybx_bnd, current%fyby_bnd, current%fybz_bnd, current%fyen_bnd, current%fypo_bnd)

!         case (4) ! to be continued
!            call mlw (current%ro, current%mx, current%my, current%mz,                                                              &
!               current%bx, current%by, current%bz, current%en, current%az, current%po, current%gx, current%gy,                     &
!               nx_block, ny_block, ng, n_neighs, current%dx, current%dy, min_dt, gm, cd, cfl, div_clean_flag,                      &
!               current%fxro_bnd, current%fxmx_bnd, current%fxmy_bnd, current%fxmz_bnd,                                             &
!               current%fxbx_bnd, current%fxby_bnd, current%fxbz_bnd, current%fxen_bnd, current%fxpo_bnd,                           &
!               current%fyro_bnd, current%fymx_bnd, current%fymy_bnd, current%fymz_bnd,                                             &
!               current%fybx_bnd, current%fyby_bnd, current%fybz_bnd, current%fyen_bnd, current%fypo_bnd)

   !      case (4)
   !         call ucs (current%ro, current%mx, current%my, current%mz,                                                                 &
   !            current%bx, current%by, current%bz, current%en, current%gx, current%gy,                                                &
   !            nx_block, ny_block, ng, n_neighs, current%dx, current%dy, min_dt, gm,                                                  &
   !            current%fxro_bnd, current%fxmx_bnd, current%fxmy_bnd, current%fxmz_bnd,                                                &
   !            current%fxbx_bnd, current%fxby_bnd, current%fxbz_bnd, current%fxen_bnd, current%fxpo_bnd,                              &
   !            current%fyro_bnd, current%fymx_bnd, current%fymy_bnd, current%fymz_bnd,                                                &
   !            current%fybx_bnd, current%fyby_bnd, current%fybz_bnd, current%fyen_bnd, current%fypo_bnd)
         end select

         if (ns .eq. 1 .and. gravity_reset_flag .eq. 1) then
            do j = 1 + ng, ny_no_bnd + ng
            do i = 1 + ng, nx_no_bnd + ng
               current%gx(i, j) = current%gx(i, j) - current%mx(i, j) / current%ro(i, j) / min_dt
               current%gy(i, j) = current%gy(i, j) - current%my(i, j) / current%ro(i, j) / min_dt
            enddo
            enddo
!if (cpu_id .eq. 0 .and. all(reshape(current%id,(/2/)) .eq. (/1,1/))) then
!print*, current%id
!print*, current%gy(16,:)
!endif

            do j = 1, ny_block
            do i = 1, nx_block
               current%ro(i, j) = ro_old(i, j)
               current%mx(i, j) = mx_old(i, j)
               current%my(i, j) = my_old(i, j)
               current%mz(i, j) = mz_old(i, j)
               current%bx(i, j) = bx_old(i, j)
               current%by(i, j) = by_old(i, j)
               current%bz(i, j) = bz_old(i, j)
               current%en(i, j) = en_old(i, j)
               current%az(i, j) = az_old(i, j)
               current%po(i, j) = po_old(i, j)
            enddo
            enddo
         endif
      enddo

      if (damping_flag .eq. 1) then
         call damping_zone(current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en,         &
            current%ro0, current%mx0, current%my0, current%mz0, current%bx0, current%by0, current%bz0, current%en0,                &
            current%x, current%y, nx_block, ny_block)
      endif

!      if (div_clean_flag .eq. 2) then
!         call mpact (et, current%ro, current%mx, current%my, current%bx, current%by, current%en, current%po,                       &
!            ro_old, mx_old, my_old, current%dx, current%dy, min_dt, nx_block, ny_block)
!      endif

   endif

   parent => current%parent
   if (parent%lv .ge. 1 .and. parent%cpu_id .ne. cpu_id) then
      parent_id_p1 = parent%cpu_id + 1
      n_blocks_update_send(parent_id_p1) = n_blocks_update_send(parent_id_p1) + 1
      blocks_update_info_send(:, :, n_blocks_update_send(parent_id_p1), parent_id_p1) = current%id
      if (.not. parent%var_allocate_flag) then
         allocate(parent%ro(nx_block, ny_block), parent%mx(nx_block, ny_block), parent%my(nx_block, ny_block),                     &
                  parent%mz(nx_block, ny_block), parent%bx(nx_block, ny_block), parent%by(nx_block, ny_block),                     &
                  parent%bz(nx_block, ny_block), parent%en(nx_block, ny_block), parent%gx(nx_block, ny_block),                     &
                  parent%gy(nx_block, ny_block), parent%po(nx_block, ny_block), parent%az(nx_block, ny_block))
         allocate(parent%fxro_bnd(ny_block, n_neighs), parent%fxmx_bnd(ny_block, n_neighs),                                        &
                  parent%fxmy_bnd(ny_block, n_neighs), parent%fxmz_bnd(ny_block, n_neighs),                                        &
                  parent%fxbx_bnd(ny_block, n_neighs), parent%fxby_bnd(ny_block, n_neighs),                                        &
                  parent%fxbz_bnd(ny_block, n_neighs), parent%fxen_bnd(ny_block, n_neighs),                                        &
                  parent%fyro_bnd(nx_block, n_neighs), parent%fymx_bnd(nx_block, n_neighs),                                        &
                  parent%fymy_bnd(nx_block, n_neighs), parent%fymz_bnd(nx_block, n_neighs),                                        &
                  parent%fybx_bnd(nx_block, n_neighs), parent%fyby_bnd(nx_block, n_neighs),                                        &
                  parent%fybz_bnd(nx_block, n_neighs), parent%fyen_bnd(nx_block, n_neighs),                                        &
                  parent%fxpo_bnd(ny_block, n_neighs), parent%fypo_bnd(nx_block, n_neighs))
         if (ini_value_save_flag .ne. 0) then
            allocate (parent%ro0(nx_block, ny_block), parent%mx0(nx_block, ny_block), parent%my0(nx_block, ny_block),              &
               parent%mz0(nx_block, ny_block), parent%bx0(nx_block, ny_block), parent%by0(nx_block, ny_block),                     &
               parent%bz0(nx_block, ny_block), parent%en0(nx_block, ny_block), parent%az0(nx_block, ny_block))
         endif
         parent%var_allocate_flag = .true.
      endif
   endif

   current%dt = min_dt
   current%t = current%t + min_dt

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine integrate
